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Abstract 

The  purpose  of  this  study  is  to  find  optimal  low 

thrust  capture  trajectories  in  the  restricted  three  body 

problem.  The  region  of  phase  space  which  corresponds  to 

capture  by  a  body  is  bounded  by  the  curve  of  zero  velocity 

n 

passing  through  the  equilibrium  point  (Lagrange  point)  L2. 

n 

If  the  spacecraft  is  driven  to  rest  at  L2,  capture  has  been 
achieved  by  adding  the  least  amount  of  energy.  An  optimal 
control  law  is  developed  to  achieve  this  based  on 
maximizing  the  Jacobi  integral.  The  dynamics  are  then 
linearized  around  L2  for  the  Earth^Moon  system  and  the 
shooting  method  is  used  to  solve  the  two  point  boundary 
value  problem.  This  problem  was  found  to  be  singular  so 
the  shooting  method  was  modified  to  avoid  making 
corrections  along  the  null  eigenvector.  Four  trajectories 
were  found  which  demonstrate  the  optimal  control  law 
successfully  caused  the  spacecraft  to  be  captured  by  the 
Moon.  'fl^S *  f\\i~ YUa 
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OPTIMAL  LOW  THRUST  TRAJECTORIES  FOR  PLANETARY  CAPTURE 


JU  lnt-EQflu.C-ti.Qn 

Optimal  trajectories  for  low  thrust  spacecraft  have 
been  the  subject  of  study  since  the  early  1960 's.  The 
majority  of  the  previous  studies  concentrated  on  finding 
optimal  transfer  orbits  in  the  classical  two  body  problem. 
This  study  will  be  concerned  with  locating  optimal  capture 
trajectories  for  low  thrust  spacecraft  in  the  restricted 
three  body  problem.  For  this  study,  a  low  thrust 
spacecraft  is  defined  as  having  engines  which  produce  less 
thrust  than  the  force  of  local  gravity. 

When  the  U.S.  space  station  becomes  operational,  the 
use  of  low  thrust  vehicles  will  become  a  practical  means  of 
interplanetary  travel.  Ion  engines,  which  produce  constant 
low  thrust  at  high  specific  impulse,  are  ideal  for 
transporting  large  payloads  between  planets. 

The  problem  of  interplanetary  travel  may  be  viewed  as 
a  series  of  capture  problems.  The  purpose  of  this  study  is 
to  find  optimal  low  thrust  capture  trajectories  in  the 
restricted  three  body  problem.  The  equations  of  motion  for 
a  thrusting  spacecraft  in  the  restricted  three  body  problem 
are  developed  and  the  unforced  motion  is  examined  to 
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are  developed  and  the  unforced  notion  is  examined  to 
determine  a  criterion  for  capture. 

For  the  unforced  motion,  an  integral  of  the  motion, 
the  Jacobi  integral  is  introduced.  From  this  integral  of 
the  motion,  curves  of  zero  velocity  are  defined.  These 
curves  determine  if  motion  is  possible  given  the  value  of 
the  Jacobi  integral.  The  region  of  phase  space  which 
corresponds  to  capture  by  a  body  is  bounded  by  the  curve  of 
zero  velocity  passing  through  the  equilibrium  point 
(Lagrange  point)  L2.  This  is  the  point  where  the  two 
capture  regions  touch.  If  the  spacecraft  is  driven  to  rest 
at  L2 ,  capture  has  been  achieved  by  adding  the  least  amount 
of  energy. 

An  optimal  control  law  is  developed  to  achieve  this 
based  on  maximizing  the  Jacobi  integral.  The  resulting  two 
point  boundary  value  problem  has  no  closed  form  solution, 
so  a  numerical  solution  is  developed.  The  dynamics  are 
then  linearized  around  L2  for  the  Earth~Moon  system  and  the 
shooting  method  is  used  to  solve  the  two  point  boundary 
value  problem. 

This  problem  was  found  to  be  singular  so  the  shooting 
method  was  modified  to  avoid  making  corrections  along  the 
null  eigenvector.  Four  cases  are  investigated.  Case  1 
demonstrated  the  modified  shooting  method  was  able  to  solve 
the  two  point  boundary  value  problem  if  the  initial 
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conditions  were  very  close  to  the  nominal  initial 
conditions.  Case  2  proved  the  problem  was  indeed  singular. 

Case  3  drove  a  spacecraft  which  was  spiralling  out 
from  Earth  to  rest  at  L2 .  Case  4  showed  that  the 
linearized  dynamics  were  valid  for  Case  3.  Therefore,  the 
optimal  control  law  successfully  caused  the  spacecraft  to 
be  captured  by  the  Moon. 


Theory 


In  order  for  a  spacecraft  to  be  captured  by  a 
celestial  body,  it  must  enter  the  region  of  phase  space 
where  the  gravity  field  of  the  target  body  is  dominant.  In 
the  simplest  case,  there  are  two  celestial  bodies 
(primaries)  and  the  spacecraft;  the  three  body  problem. 

With  one  celestial  body  and  a  spacecraft,  the  celestial 
body's  gravity  field  dominates  all  of  the  phase  space. 
Hence,  the  three  body  problem  is  the  simplest  nontrivial 
capture  problem.  Fortunately,  the  restricted  three  body 
problem  adequately  describes  most  capture  scenarios  in  our 
solar  system. 

The  solution  to  the  problem  of  how  to  optimally  cause 
a  low  thrust  vehicle  to  be  captured  involves  developing 
equations  of  motion,  analysis  of  the  motion,  determination 
of  the  optimal  control  law,  and  solution  of  the  two  point 
boundary  value  problem. 

The  Restricted  Three  Bodv_ Problem 

This  study  will  be  confined  to  the  restricted  three 
body  problem.  The  restricted  three  body  problem  consists 
of  two  bodies  revolving  about  their  center  of  mass  in  a 
circular  orbit  at  a  distance  1  from  each  other.  A  third 
body  of  negligible  mass,  moves  in  the  plane  defined  by  the 
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The  position  vector  of  the  third  body  expressed  in  the 
rotating  frame  is 

r  =  xt  +  y}  (l) 

where  1  is  a  unit  vector  in  the  x  direction 
}  is  a  unit  vector  in  the  y  direction. 

The  velocity  is  found  by  applying  Coriolis'  theorem  for 
taking  derivatives  in  a  rotating  frame,  which  is 

dr  dr 

v  =  -  =  -  +  u  x  r  (2) 

dt  i  dt  R 

where 

dr 

-  is  the  inertial  time  derivative 

dt  z 

dr 

-  is  the  time  derivative  in  the  rotating  frame 

dt  R 

u>  is  the  angular  velocity  of  the  rotating  frame. 
Performing  the  cross  product  produces 

v=(x-ny)t+(y+nx)^  (3) 

Therefore,  the  kinetic  energy  of  the  third  body  is 

T  =  ^m(vv)  =  *5m[(*  -  fly)2  +  (y  +  nx)2]  (4) 

The  mass  of  the  third  body  plays  no  role  in  the  motion, 
hence  the  kinetic  energy  per  unit  mass  becomes 

T/m  =  ^[(x  -  fly)2  +  (£  +  nx)2]  (5) 

The  coordinates  x  and  y  may  be  nondimensionalized  such  that 
the  new  coordinates  x  and  y  equal  the  old  x  and  y  divided 
by  the  distance  between  the  primaries,  1  =  a  +  b,  hence 

T/ml2  =  % [ ( x  -  ny)2  +  (y  +  Ox)2]  (6) 
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Time  is  nondimens ionali  zed  by  dt*  =  ndt,  where  t*  is  the 

nondimensional  time.  Therefore, 

dx  dx 

-  =  n  - 

dt  dt* 

Denoting  x  =  dx/dt*  and  y  =  dy/dt*,  produces 

T/ml2  =  ^fi2(x2  -  2ky  +  y2  +  y2  +  2x£  +  x2 ) 

Kinetic  energy  is  nondimensionalized  by  dividing  by  n2  to 
yield 

T*  =  T/ml2n2  =  ^ { ( x2  +  $2)  +  2(xy  -  xy)  +  (x2  +  y2)}  (7) 
The  potential  energy  of  the  third  body  is  the  sum  of 
the  classical  gravitational  potentials  due  to  the  two 
primaries,  so  the  potential  energy  expressed  in  dimensional 
coordinates  is 

V  =  -Gmfmj/r!  +  m2/r2)  (8) 

where 

G  is  the  universal  gravity  constant 

m-^  and  m2  are  the  masses  of  the  primaries 


and 


rx  =  [(x  -  b)2  +  y2]* 

(9) 

r2  -  t(x  +  a)2  +  y2]1* 

(10) 

The  potential 

energy  per  unit  mass  is  given  by 

V/m  =  -Gfmi/r!  +  m2/r2) 

(11) 

Balance  between  gravitational  and  centrifugal  forces 

requires : 

Gmim2/12  =  m2afl2  =  m^bn2 

(12) 
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hence 


Gm^  =  an2l2 

(13) 

Gm2  =  bn2l2 

(14) 

G(m1  +  m2)  =  n2l2 

(15) 

Two  mass  parameters,  m  and  |i2  are  defined  as 

»1 

m  = 

m^  +  m2 

m2 

^2  = 

®1  ®2 

follows: 

and  since  the  coordinate  frame  is  centered  at 

mass  of  m^  and  m2 ,  1  and  n2  become 

the  center  of 

m  =  a/1 

(16) 

\i2  =  b/1 

Substituting  (13)  and  (14)  into  (11)  leads  to 

(17) 

V/m  =  -(an2l2/ri  +  bn2l2/r2) 
Dividing  by  n2  yields 

(18) 

V/mfl2  =  -(al2/ri  +  bl2/r2) 

(19) 

Now  r!  and  r2  are  nondiraensionalized  by  dividing  by  1. 

Then  substituting  (16)  and  (17)  into  (9)  and  (10),  the 
nondiinensional  r*  and  r2  become 

ri  *  [(x  -  n2)2  +  y2]3*  (20) 

r2  =  [(x  +  m)2  +  y2]^  (21) 

and  therefore  (19)  becomes 

V/ran2  =  -(al/r^  +  bl/r2) 


Substituting  for  a  and  b  from  (16)  and  (17)  and  dividing  by 
l2  yields  the  nondimensional  form  of  the  potential: 

v*  =  v/mn2l2  =  -(m/ri  +  m/r2)  (22) 

Hamilton's  eguations  are  now  used  to  derive  four  first 
order  equations  of  motion  from  the  expressions  for  kinetic 
and  potential  energy.  Two  new  coordinates  are  introduced, 

m 

called  conjugate  momenta,  denoted  Px  and  py.  Since  the 
potential  is  not  explicitly  dependent  on  the  rates  x  and  y, 
the  conjugate  momenta  may  be  found  by 

px  =  3T*/0x 

Py  «  3T*/dy 

Taking  the  partials,  px  and  Py  become 

px  =  x  -  y  (23) 

Py  =  y  +  x  (24) 

The  kinetic  energy  is  separated  into  terms  which 
involve  the  rates,  x  and  y,  quadratrically ,  linearly,  and 
those  terms  which  do  not  involve  the  rates;  T*2,  T*lr  and 
T*q  respectively.  The  Hamiltonian  is 

H  =  T*2  +  V*  -  T*0 

which  yields 

H  =  *(x2  +  $2)  -  (\i1/r1  +  n2/r2)  -  *(x2  +  y2)  (25) 

Substituting  for  x  and  y  using  (23)  and  (24)  leads  to 

H  ■  ^(Px2  +  2PXy  +  Py2  "  2Pyx)  “  (HiAi  +  H2/r2)  (26) 


If  the  only  force  acting  is  gravity,  Hamilton's 
equations  for  this  system  are 

x  =  3H/3px  }  =  3H/3py  (27) 

px  =  -3H/ax  Py  =  -aH/ay 

The  thrust  vector  applied  to  the  spacecraft  is  defined  as 

/ 

follows : 


P  =  F  cos S  i  +  F  sins  3  (28) 

where  F  is  a  nondimensional  constant  thrust  magnitude,  (F  = 
thrust/min2 ) ,  and  s  is  the  angle  of  the  thrust  vector 
relative  to  the  local  horizontal.  This  is  a  force  which  is 
not  accounted  for  in  the  potential,  so  Hamilton's  equations 
become 

X  =  3H/apx  y  =  3H/3py  (29) 

£x  =  -3H/3X  +  Qx  Py  =  -3H/3y  +  Qy 

where 


qx  =  p»av/ax  Qy  =  p»av/ay 

which  yields 

Qx  =  F  cos S  Qy  =  F  sins  (30) 

Substituting  (30)  into  Hamilton's  equations  and  taking  the 
indicated  derivatives  of  H,  the  equations  of  motion  become 

x  =  px  +  y  (31) 

y  -  Py  -  x 

Px  =  Py  "  tUi(x  -  U2)/rl3  +  M-2 ( x  +  ^l)/r23]  +  F  cosS 
£y  *  "Px  “  y[Hi/ri3  +  H2/r23]  +  F  sins 
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Note  that  iij.  and  n2  maY  be  between  0  and  1  and  112  =  1 
-  nj_,  so  these  two  parameters  may  be  consolidated. 

Let  11  *  n2>  then  (i^  =  1  -  |j. ,  which  leads  to  the  equations 
of  motion: 

x  =  px  +  y  (32) 

*  *  Py  "  x 

Px  =  Py  ~  C(l-tO(x  -  n)/rx3  +  ji(x  +  (l-ji))/r23]  +  F  cos# 
Py  =  "Px  ‘  y[(l“>i)/ri3  +  n/r23 ]  +  F  sin# 

where 

r1  =  [(x  -  n)2  +  y2]^ 
r2  =  [(x  +  (1-n))2  +  y2]^ 

Regions  of  Motion  and  Capture  Criterion 

The  equations  of  motion  for  a  thrusting  spacecraft  in 
the  restricted  three  body  problem  are  shown  in  eqn.  (32). 
Unforced  motion  in  the  restricted  three  body  problem  has 
been  studied  since  1772.  (4:4)  Analysis  of  this  motion 
will  help  determine  a  definition  of  capture  and  a 
performance  index  for  an  optimal  control  law. 

For  the  unforced  motion,  all  forces  are  accounted  for 
in  the  Lagrangian.  Since  the  Hamiltonian  is  not  an 
explicit  function  of  time,  ft  =  3H/3t  =0,  it  is  an  integral 
of  the  motion,  called  the  Jacobi  integral.  The 
introduction  of  this  integral  was  Jacobi's  major 
contribution  to  this  problem  in  1836.  (4:4)  Since  ft  »  0,  H 
is  constant,  therefore  the  motion  must  take  place  in  such  a 
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manner  which  renders  the  Hamiltonian  constant  over  time. 
This  leads  to  some  deductions  about  the  motion. 

Recall 

H  =  +  y2)  -  (UiAx  +  n2/r2)  "  *(x2  +  y2) 

which  is  constant  over  time.  Therefore,  setting  H  =  -  C  (a 
constant)  and  rearranging  yields 

x2  +  y2  =  2(^1/r1  +  ji2/r2)  +  (x2  +  y2)  “  2C 
which  may  be  written 

X2  +  y2  *  2U(x,y )  -  2C  (33) 

The  left  hand  side  of  (33)  is  always  greater  than  zero, 
which  implies  motion  only  exists  if  U(x,y)  >  C.  The 
equation  U(x,y)  =  C  defines  "curves  of  zero  velocity." 

Figs.  2  through  5  depict  these  curves  for  various  \i.  The 
constant  C  is  determined  by  the  initial  conditions. 

Setting  y  *  0  determines  how  the  Hamiltonian  or  Jacobi 
integral  varies  along  the  x-axis  for  various  values  of  ^ 
(shown  in  Fig.  6).  There  are  three  local  maxima  along  the 
x-axis.  These  are  unstable  equilibrium  points  called  L^, 

L2 ,  and  L3,  which  are  referred  to  as  Lagrange  points. 
(4:142) 

Curves  of  zero  velocity  provide  information  about  the 
possible  motion  in  regions  of  the  x-y  plane.  In  order  to  be 
captured  by  a  celestial  body,  a  spacecraft  must  enter  the 
region  of  phase  space  where  the  motion  consists  of  a  closed 
orbit  about  that  body.  It  is  difficult  to  visualize  a  four 


12 


dimensional  phase  space,  however,  the  curves  of  zero 
velocity  characterize  the  motion  in  terms  of  the  two 
position  coordinates  and  the  Jacobi  integral  (which  is  a 
measure  of  energy). 

The  region  of  the  phase  space  which  defines  capture  is 
bounded  by  the  contour  of  the  2ero  velocity  curve  passing 
through  L2  because  if  the  spacecraft  is  in  this  region  of 
the  x-y  plane  with  a  value  of  the  Jacobi  integral  less  than 
that  of  the  contour  passing  through  L2,  the  spacecraft 
cannot  leave  this  region.  This  is  the  definition  of 
capture . 

The  question  now  becomes:  what  direction  should  the 
thrust  be  applied  to  optimally  cause  the  spacecraft  which 
is  captured  by  one  body  to  be  captured  by  the  other?  The 
two  capture  regions  touch  at  L2.  L2  is  a  saddle  point 
because  the  Jacobi  integral  has  a  local  maximum  at  L2  along 
the  x  direction,  but  has  a  local  minimum  along  the  y 
direction  there.  To  enter  the  target  capture  region  at  any 
point  other  than  L2  requires  a  Jacobi  integral  (and  hence 
energy)  greater  than  that  of  the  zero  velocity  curve 
passing  through  L2  because  the  spacecraft  must  leave  the 
initial  capture  region.  Therefore,  this  is  the  point  where 
a  spacecraft  captured  by  one  body  may  be  caused  to  be 
captured  by  the  other  by  adding  the  least  energy. 
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Fig.  2.  Curves  of  Zero  Velocity  for  ix 
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The  general  optimal  control  problem  is  to  determine 
the  control  which  will  drive  a  system  to  a  desired  set  of 
final  conditions,  while  minimizing  a  cost  function  or 
performance  index.  This  cost  function  is  usually  expressed 
in  the  form  of  an  integral 

PI  =  L(x, u, t )  dt 

•  ti 

where  x  is  the  state  vector  and  u  is  the  control  vector. 

The  Jacobi  integral  suggests  itself  as  a  logical 
choice  for  a  performance  index  because  the  Jacobi  integral 
must  be  increased  to  cause  capture.  A  large  difficulty 
with  optimal  control  theory  is  choosing  a  proper 
performance  index  since  a  system  may  behave  optimally  for 
the  given  performance  index,  there  is  no  guarantee  the 
given  performance  index  is  really  "optimal."  Therefore, 
performance  indices  with  physical  significance  usually  are 
preferred.  In  this  case,  the  Jacobi  integral  or 
Hamiltonian  has  a  definite  physical  significance. 

The  control  variable  for  this  problem  is  9  .  Constant 
thrust  will  be  assumed  to  be  applied  to  a  spacecraft  of 
constant  mass.  Since  the  thrust  magnitude  is  constant,  the 
problem  is  to  determine  the  function  0(t)  which  will  drive 
the  system  to  be  captured  as  quickly  as  possible. 
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The  equations  of  motion  have  been  developed,  a 
performance  index  chosen,  and  desired  final  conditions 
determined.  The  calculus  of  variations  is  used  to 
determine  necessary  conditions  to  minimize  the  performance 
index.  In  order  to  drive  the  spacecraft  toward  L2  from  the 
vicinity  of  one  of  the  bodies,  the  Jacobi  integral  must  be 
increased.  Thus  if  -H  is  minimized,  increase  in  H  will  be 
maximized,  driving  the  spacecraft  toward  capture.  This 
performance  index  expressed  in  integral  form  becomes 

r*. 

PI  =  -H  dt 

Jti 

Since  the  system  is  constrained  to  obey  the  equations  of 
motion,  Lagrange  multipliers  (costates)  are  introduced  so 
the  performance  index  becomes 

'U 

PI  -  (-H  +  ATf )dt 

•  tj 

where  A  is  a  vector  of  costates  and  the  elements  of  f  are 
the  equations  of  motion  in  vector  form. 

The  chain  rule  is  used  to  determine  H,  which  is 
H  *  x  3H/3X  y  3H/3y  +  £x  3H/apx  +  Py  3H/3Py 
The  partial  derivatives  are  found  by  using  Hamilton's 
equations  (29)  so  that  H  becomes 

H  =  (F  cos*  -  Px)x  +  (F  sin*  -  Py)y  +  xpx  +  ypy 
which  simplifies  to 

H  *  £  Fcos*  +  y  Fsin*  (34) 
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The  necessary  conditions  for  a  minimum  are  found  by 
taking  the  first  variation  of  PI,  integrating  by  parts,  and 
setting  that  equal  to  zero,  which  produces  the  Euler 
equation  and  boundary  conditions.  This  is  presented  by 
Schultz  and  Melsa  for  a  general  PI.  (3:233-249)  Those 
results  will  be  employed  to  determine  the  optimal  control 
law.  The  state  function  of  Pontryagin  (also  known  as  the 
control  Hamiltonian)  is 

Hc  =  -H  +  ATf (x,9 ,t) 
where  the  equations  of  motion  are 

x  =  f  (x,9 ,t) 


and 

x  =  (x  y  px  py)T 

and  A  is  a  vector  of  four  Lagrange  multipliers  or  costates. 
Substituting  for  H  and  f,  the  control  Hamiltonian  becomes 
Hc  *  “XF  cos*  -  yF  sin*  +  AjJp*  +  y)  +  A2(py  -  X)  + 
*3<Py  “  t(l-U)(x  -  u)/rx3  +  n(x  +  (l-|i))/r23]  +  F  cos*} 

+  a4<“Px  _  YtU-fcO/ri3  +  H/r23  ]  +  F  sin*  }  (35) 
The  optimal  control  is  found  by  setting  3 Hc/3*  =  o,  which 
yields 

3 Hc/3  9  =  xF  sin*  -  yF  cos*  -  A3F  sin*  +  a4F  cos*  =  0 
Dividing  through  by  F  and  grouping  sin*  and  cos*  leads  to 
(x  -  >3)sin*  +  (a4  -  y)cos*  =  0 
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Hence,  the  optimal  control  law  is 

(y  -  *4) 

tan#  =  - 

(x  -  x 3 ) 

and  therefore  the  optimal  thrust  angle  $*  is 
8*  =  tan_1[(y  -  x4)/(x  -  x3)] 

The  sine  and  cosine  of  the  optimal  thrust  angle  are 

y  -  x4 


sm  e  = 


cos 


[(y  -  A4)2  +  (*  - 

X  -  X3 

[(y  -  * 4 ) 2  +  (x  - 

Substituting  for  x  and  y  from  (32)  produces 

Py  -  x  -  x4 


sin  6*  = 


cos  8  = 


( ( Py  -  x  -  a4)2  +  (px  +  y  -  x 3 ) 2 ] ^ 
Px  +  y  -  a3 


(36) 


(37) 


(38) 


(39) 


((py  -  X  -  X4)2  +  (px  +  y  -  ^  3 ) 2 ] ^ 

The  optimal  trajectory  will  obey  the  following: 

X  =  3HC/3A  I  =  -3HC/3X 

The  first  four  equations,  x  =  3HC/3X,  are  the  equations  of 
motion  (32)  with  the  optimal  thrust  angle  8*  substituted 
for  6.  The  remaining  four  equations,  which  describe  how 
the  costates  vary  in  time,  are  as  follows: 


lx  -  -F  sin 8*  +  x: 


x3(l  -  n)  HA  3 
+  -  +  - 

rl3  r23 

3(1  -  H ) ( x  -  H)[A3(x  -  H)  +  A4y] 
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(40) 


3u[x  +  (l  -  n)](x3[x  +  (l  -  H)]  +  x4y) 


r2' 


.  *4<1  -  H)  U*4 

=  F  COS?  -  Ai  +  -  +  - 


ri  r2 

3y(l  -  »4)[x3(x  -  ia)  +  x4y] 


3ny{x3[x  +  (1  -  p.)  ]  +  x4y) 


r2 


=  F  cos?  -  +  A, 


(41) 

(42) 

(43) 


x4  =  F  sin?*  -  x 2  -  x 3 
The  second  derivative  of  Hc  with  respect  to  ?  must  be 
greater  than  0  to  insure  a  local  minimum.  (2:184). 
Therefore,  for  a  local  minimum 

a2 Hc/a?2  =  (x  -  x3)cos?  -  (x4  -  y)sin?  >  o  (44) 
Along  a  stationary  trajectory,  cos?  and  sin?  are  given  by 
(38),  which  leads  to 


(x  -  x3)2  +  (y  -  a4)2 


[(py  -  X  -  a4)z  +  (px  +  y  -  x3)<] 


2 


>  0 


(45) 


which  is  satisfied  because  the  positive  sguare  root  is 
chosen  by  the  definition  of  the  sine  and  cosine,  so  along 
any  stationary  path,  a  local  minimum  is  guaranteed. 
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Since  the  equations  of  motion  and  optimal  control  law 
have  been  developed,  solving  the  boundary  value  problem  is 
next  step.  Unfortunately,  boundary  conditions  are 
determined  at  both  the  initial  and  final  times.  At  the 
initial  time,  t^,  the  four  states  are  known  but  the 
costates  are  not.  Also,  at  the  final  time,  tf,  the  desired 
states  are  determined  and  the  costates  are  undetermined. 

Taking  the  first  variation  of  the  performance  index, 
PI,  and  integrating  by  parts  from  t^  to  tf  yields  the 
following  boundary  condition: 

-Aj^dx  -  A2dy  -  A3dpx  -  A4dpy  +  Hcdt  =  0  (46) 

evaluated  at  tf.  Since  the  states  at  tf  are  fixed, 

dx  *  dy  =  dpx  =  dPy  =  o 

so  Hcdt  =0.  If  tf  is  specified,  dt  =  0  so  the  boundary 
conditions  are  satisfied.  If  tf  is  not  specified,  tf  is 
chosen  so  that  Hc  =  0  at  tf.  This  value  of  tf  minimizes 
PI. 

The  problem  is  to  determine  the  initial  costates  to 
drive  the  system  to  rest  at  L2.  No  closed  form  solution  to 
this  problem  is  known  since  no  closed  form  solution  to  the 
restricted  three  body  problem  is  known  and  a  nonlinear 
control  is  being  applied,  therefore  a  numerical  solution  is 
called  for.  The  shooting  method  was  chosen  to  solve  this 
problem. 


The  shooting  method  is  derived  from  the  behavior  of 
linearized  dynamics.  The  equations  of  motion  form  a  set  of 
eight  nonlinear  ordinary  differential  equations  which  is 
expressed  in  vector  form  as 

X  =  g(X,t)  (47) 

Expanding  about  a  desired  trajectory,  Xq,  in  a  Taylor 
series  by  letting  X  =  Xq  +  6X  and  substituting  into  the 
equations  of  motion  produces 
«  d(XG  +  SX) 

X  =  -  =  g(X0  +  6X,t)  (48) 

dt 

dXQ  d4X 

-  +  -  =  g(X0,t)  +  Vxg(X0,t)6X  +  0(2)  (49) 

dt  dt 

Cancelling  the  equations  of  motion  for  the  desired 
trajectory,  yields  the  equations  of  variation: 

SX  *  A( t ) £X  (50) 

where  A(t)  =  Vxg(X0,t),  which  is  evaluated  along  the  time 
dependent  desired  trajectory.  The  equations  of  variation 
are  linear  ordinary  differential  equations,  so  a  set  of 
independent  solutions  may  be  constructed  as  follows: 

^i  “  A(t)^i  ?  ^ij(t^)  =  (51) 

where  is  the  Kroenecker  delta.  Any  solution  to  the 
equations  of  variation  may  be  written 

SX(t)  =  2  4Xi(ti)^i(t) 

where  4Xi(tjJ  are  the  components  of  4X  at  t^.  The  form 
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a  fundamental  set  of  solutions,  so  assembling  the  vector 
columns  of  <t>±  into  a  matrix  *  leads  to  the  expression: 

6X(t)  =  *(t,ti)iX(ti)  (52) 

The  state  transition  matrix,  •  satisfies 

*  =  A(t)*  ;  #(ti,ti)  =  I  (53) 

where  I  is  the  identity  matrix.  This  equation,  along  with 
the  equations  of  motion  are  numerically  integrated  from  tj_ 
to  tf.  (5:24-27) 

Applying  this  to  the  case  at  hand,  X  is  defined  as 


Eqn.  (52)  is  partitioned  as  follows: 

|fix(tf)|  ^  T  *xx  |  *XA  6x(t^)|  ^ 

UMtf)J  l  #AX  I  *AA  JUx(ti)J 

Hence,  the  variation  in  the  final  states  due  to  small 
changes  in  initial  conditions  may  be  written: 

6x(tf)  =  *xx6x(ti)  +  *xA6A(ti)  (55) 

but  since  x(tj^)  is  fixed,  4x(t^)  =  (0)  so 

Sx(tf)  =  *xASA(ti)  (56) 

Eqn.  (56)  is  an  expression  which  determines  how  to  change 
A ( tj[ )  to  obtain  a  desired  change  in  x(tf).  Solving  for 
4a  ( t j_ ) ,  produces 

4A(ti)  =  #XA*1«X(tf)  (57) 

The  shooting  method  starts  by  estimating  a  set  of 
initial  costates  and  integrating  the  equations  of  motion 


and  the  equation  for  *,  (53),  from  t^  to  tf.  This  produces 
x(tf)  for  the  first  iteration,  £x(tf)  is  the  difference 
between  x(tf)  and  the  desired  x(tf).  Then  equation  (57) 
determines  the  correction  to  the  initial  costates.  The 
next  set  of  initial  costates  will  be  the  previous  initial 
costates  minus  6A(t^),  and  the  process  is  repeated  until  a 
trajectory  from  x(tj^)  to  the  desired  x(tf)  is  found. 
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Ill, 


The  first  practical  application  of  this  theory  will  be 
in  the  Earth-Moon  system.  This  case  has  the  highest  value 
of  pl  for  any  three  body  problem  in  the  solar  system.  At 
this  point  attention  is  confined  to  the  Earth-Moon  case  for 
fixed  final  times.  For  the  Earth-Moon  system: 

n  =  2.649  X  10"6  radians/sec  (1:337) 

1  *  384,400  km.  (1:323) 

4  -  0.0121506683 

and  L2  is  located  at  x^=  -0.8369147190  along  the  x-axis. 
(4:221) 

The  first  step  was  to  verify  the  equations  of  motion 
and  the  numerical  integrator.  The  integration  algorithm 
used  was  a  routine  called  naming,  which  is  used  in  Mech. 
731.  (5)  The  equations  of  motion  (eqns. 

(32) ,(38) ,(40) -(43))  were  programmed  and  then  tested. 

There  is  a  periodic  orbit  which  alternates  between  the 
earth  and  the  moon  found  by  Davidson  in  1964  for  4  = 

1/81.45  and  no  thrust.  (4:509-510)  This  orbit  was  computed 
with  initial  conditions:  x  =  -1.15,  y  *  0.0,  px  =  0.0,  py 
»  -1.15868829;  with  a  step  size  of  1000  integration  steps 
per  time  unit.  This  value  of  step  size  was  chosen  to 
insure  that  H  remained  constant.  The  orbit  was 
successfully  reproduced,  (shown  in  Fig.  7),  verifying  the 
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Evaluating  at  L2;  x  =  x^  y  =  0,  £  =  0,  y  =  0,  leads  to 


A(t)  = 


1  0 


0  1 


0  1 


0  a42  -1  0 


where 


a31  =  10.29519541 
a42  =  -5.147597707 

A  new  set  of  variables  are  defined  as  follows: 


*  '  *La 


l  x4  J  l  Py  -  xLjlJ 

Therefore,  the  equations  of  motion  near  L2  with  thrust 
applied  become 

Xi  =  x2  +  x3  (61) 

x2  »  -X3  +  x4 
x3  =  831X1  +  x4  +  FcosS 
x4  *  a42x2  -  x3  +  Fsintf 

Using  Hamilton's  equations  and  the  chain  rule  as  in  (34),  H 
becomes 

H  =  xiFcos^  +  x2Fsine  (62) 

Substituting  for  x3  and  x2  from  (61)  leads  to 

H  =  (x2  +  x3)Fcoss  +  (-x^  +  x4)Fsin* 
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The  control  Hamiltonian  becomes 

Hc  =  -ft  +  a3(x2  +  x3)  +  A  2  ( “X1  +  x4 ) 

+  A3(a3ixi  +  x4  +  Fcosfi )  +  ^4(842X2  -  x3  +  Fsins) 
The  optimal  control  law,  a HC/<3S  *  0,  is 

{-Xi  +  x4  -  A4) 

tan?  =  -  (63) 

(x2  +  x3  -  a3) 

hence 

(-x3  +  x4  -  A4) 

sins  =  - (64) 

[(x2  +  x3  -  a3)2  +  (-x2  +  x4  -  a  4 ) 2  ]  * 

(X2  +  X3  -  a3) 

COSS  =  - (65) 

C ( X2  +  X3  -  A,)2  +  (-Xi  +  x4  -  a4)2]^ 

As  before,  the  costates  obey  A  =  -aHc/3x,  where 

A1  *  x2  ~  a31A3  "  FsinS  (66) 

A  2  =  ~A 1  -  a42A4  +  FcosS 

a3=-A2+A4+  FcosS 

a  4  *  -A 2  “  a3  +  Fsins 

The  state  transition  matrix,  *,  for  this  system  is 
found  by  (53)  where  the  matrix  A  in  this  case  is  the  8  by  8 
Jacobian  of  the  equations  (61)  and  (66).  The  eight 
equations  (61)  and  (66),  along  with  the  differential 
equation  for  ♦  may  be  integrated  from  t^  to  tf.  Then  *  may 
be  verified  by  changing  the  initial  conditions  of  the 
states  and  costates  by  a  known  amount.  The  change  in  the 
final  conditions  due  to  this  change  in  initial  conditions 
is  predicted  by  (56).  These  predicted  values  were  compared 
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to  the  actual  changes  in  final  conditions  for  several  cases 
to  verify  *  was  computed  correctly. 

Next,  the  shooting  method  was  implemented  using 
equations  using  (61)  and  (66).  However,  eqn.(57)  could  not 
be  solved  because  #xA  was  singular.  When  this  occurs, 
changes  may  be  made  to  the  initial  costates  which  cause  no 
change  in  the  final  states.  One  solution  to  this  problem 
is  to  not  make  changes  to  the  initial  costates  along  the 
eigenvector  associated  with  the  zero  eigenvalue  of  *xA . 

Recalling  (56),  6x(tf)  =  *XAtfA(tj[),  *xX  may  be 
expressed 

*xA  =  EJE"1  (67) 

where  E  is  the  matrix  of  eigenvectors  and  J  is  a  diagonal 
matrix  with  the  eigenvalues:  £]_,  £ 2>  and  ?4»  on  the 
diagonal,  so  therefore, 

«x(tf)  =  EJE~1iA(ti)  (68) 

Premultiplying  by  E”1  yields 

E_16x( tf )  =  JE“1«A(ti)  (69) 

The  following  vectors  may  be  defined 

q  =  E_16x(tf) 
s  =  E_16A(ti) 

Therefore  (69)  is  expressed 

q  =  Js  (70) 

which  is  an  uncoupled  set  of  equations: 

qi  =  ? i®i  ^  =  (71) 
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Solving  for  yields 

si  =  Qi/ £  i  (72) 

Now  if  +xA  is  singular,  at  least  one  eigenvalue  =  0. 
Since  sj  is  arbitrary,  set  Sj  =  0.  Solving  for  £>(ti) 
produces 

6\( tjj  =  Es  (73) 

Since  sj  is  set  to  zero,  no  corrections  are  made  to  the 
initial  costates  along  the  null  eigenvector. 

Case  1 

The  purpose  of  Case  1  was  to  find  any  converging 
trajectory.  The  radius  of  convergence  for  the  modified 
shooting  method  was  found  to  be  very  small.  In  order  to 
produce  a  converging  trajectory,  the  initial  states  and  the 
first  estimate  of  the  initial  costates  had  to  integrate  to 
within  about  10~3  of  L2  in  position  and  10“6  in  velocity. 
This  trajectory  was  found  by  integrating  backwards  from  1,3. 
For  Case  1,  tf  =  0.7  which  is  equal  to  73.403  hours,  and  F 
=  o.l.  The  initial  conditions  for  this  case  are  as 
follows: 

X1  *  0.015865410 
x2  =  0.018917850 
x3  *  -0.07711602 
X4  =  -0.03304073 

and  the  initial  thrust  angle  is  9  *  11.647  degrees.  The 
initial  distance  from  L2  is  9490.84  km.  At  the  initial 
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conditions,  the  magnitude  of  gravity  and  Coriolis  forces  is 
0.1318,  so  the  thrust  applied  is  less  than  gravity  and 
Coriolis.  Therefore  this  is  a  low  thrust  situation  by  the 
definition  stated  in  the  introduction.  Fig.  8  shows  x^  vs 
x2  and  Fig.  9  shows  the  thrust  angle  B  vs  time. 

The  eigenvector  associated  with  the  zero  eigenvalue  of 
*XA  for  this  trajectory  is 

i  =  (1.0  -0.12876  0.29116  0.060112)T  (74) 

The  purpose  of  Case  2  was  to  verify  that  changes  in  the 
initial  costates  along  this  vector  cause  no  change  in  the 
final  states. 

Case -2 

This  case  used  the  same  initial  states,  final  time  and 
thrust  as  Case  1.  The  initial  costates  were  changed  by  0.1 
multiplied  by  (74).  The  trajectory  for  this  case  is  shown 
in  Fig.  10  and  e  vs  time  is  shown  in  Fig.  11.  0(t)  for 

this  case  exactly  matched  0(t)  in  Case  1  so  the  same 
control  was  applied  despite  different  costates.  Since  the 
control  applied  was  the  same,  the  trajectory  was  exactly 
the  same  as  Case  1.  This  confirms  that  changing  the 
initial  costates  along  the  null  eigenvector  produced  no 
change  in  the  final  states. 
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Cass.  2 


The  purpose  of  Case  3  was  to  drive  a  spacecraft 
initially  spiralling  out  from  Earth  to  rest  at  L2.  This 
case  was  chosen  because  an  outward  spiral  is  the  optimal 
way  to  escape  the  primary  in  the  two  body  problem  for 
continuous  low  thrust.  Since  this  spacecraft  was 
spiralling  out  from  Earth,  the  initial  velocity  corresponds 
to  circular  orbit  velocity  at  the  initial  distance  from 
Earth . 

Values  for  F  and  tf  as  well  as  initial  costates  were 
refined  until  after  several  forward  and  backward 
integrations,  a  converging  trajectory  was  found.  For  this 
case,  F  =  2.4  and  tf  =  0.1044,  which  corresponds  to  10.9475 
hours.  The  initial  position  was: 

XX  =  0.9369842  X  10“3 
X2  =  0.1298685  X  10"1 

which  puts  the  spacecraft  0.848128402  units  or  326020.5577 
km.  from  Earth.  Circular  velocity  for  this  distance  is 
1.105658823  km/sec  which  is  1.085814904  units.  This  is  the 
inertial  velocity  of  the  spacecraft  at  tf  and  the  velocity 
vector  is  perpendicular  to  the  radius  vector  from  the  Earth 
to  the  spacecraft.  This  radius  vector  is  at  179.1227348 
degrees  from  the  +x  direction  at  tf,  so  the  initial 
inertial  velocity  expressed  in  the  rotating  frame  is 


-0.016624443 


Vy  =  -1.085649737 

Fortunately  the  conjugate  momenta  px  and  Py  turn  out  to  be 
the  inertial  velocities  of  the  spacecraft.  Since  x3  =  px 
and  x4  =  py  -  xL,  x3  and  x4  become 

X3  =  -0.016624443 
X4  =  -0.248728034 

Fig.  12  presents  the  trajectory  and  Fig.  13  presents  0  vs 
time  for  this  case.  Along  this  trajectory,  there  were  two 
eigenvalues  which  were  nearly  zero .  This  means  changes  to 
the  initial  costates  along  two  null  eigenvectors  produce  no 
change  in  the  final  states.  This  might  be  expected  since 
the  trajectory  is  integrated  over  a  shorter  time  span  so 
changes  in  the  control  would  have  less  effect  on  the  final 
states.  The  magnitude  of  gravity  and  Coriolis  forces  at  t^ 
is  0.2449,  (about  one  tenth  of  the  thrust  value). 

Although  this  may  not  be  low  thrust  by  strict  definition, 
it  must  be  noted  that  at  L2  all  forces  cancel,  so  the 
stated  definition  breaks  down  near  L2.  F  =  2.4  could 
actually  be  a  very  small  thrust;  for  example,  for  a 
spacecraft  of  mass  15,447  kg.,  F  =  2.4  is  equivalent  to  100 
N  of  thrust. 

Case.  4 

The  purpose  of  Case  4  was  to  demonstrate  that  the 
linearized  dynamics  were  valid  in  the  vicinity  of  L2.  For 
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Case  4,  the  initial  conditions  for  Case  3  were  translated 
back  to  the  nonlinear  states  so  that 

X  =  -0.836036896 
y  =  0.01299080 
px  =  -0.016625808 
py  «  -1.085642753 

and  the  same  initial  costates  used.  The  full  nonlinear 
dynamics  [eqns.(32),  (40)— (43) ]  were  integrated  from  these 
initial  conditions  to  tf  =  0.1044  with  F  =  2.4.  The 
trajectory  is  shown  in  Fig.  14  and  8  vs  time  is  shown  in 
Fig.  15.  Fig.  14  indicates  that  the  spacecraft  does  indeed 
go  to  1,2  as  expected.  Fig.  15  is  identical  to  Fig.  13, 
therefore,  the  control  applied  in  Case  3  and  Case  4  was 
exactly  the  same.  This  demonstrates  that  the  linearized 
dynamics  were  valid  for  Case  3. 
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Fig.  12.  Trajectory  for  Case 


Fig.  13.  Theta  (degrees)  vs  Time  for  Case 


Fig.  14.  Trajectory  for  Case 


IVjl.  Conclusions 


It  has  been  shown  that  capture  may  be  achieved  by 
driving  a  spacecraft  to  rest  at  L2  in  the  restricted  three 
body  problem.  An  optimal  control  law  was  developed  to 
achieve  this  based  on  maximizing  H.  The  resulting  two 
point  boundary  value  problem  had  no  closed  form  solution, 
so  a  numerical  solution  was  developed.  The  dynamics  were 
linearized  around  L2  for  the  Earth-Moon  system  and  the 
shooting  method  was  used  to  solve  the  two  point  boundary 
value  problem. 

This  problem  was  singular  so  the  shooting  method  was 
modified  to  avoid  making  corrections  along  the  null 
eigenvector.  Case  1  showed  that  the  modified  shooting 
method  was  able  to  solve  the  two  point  boundary  value 
problem  if  the  initial  conditions  were  very  close  to  the 
nominal  initial  conditions.  Case  2  showed  that  the  problem 
was  indeed  singular. 

Case  3  drove  a  spacecraft  which  was  spiralling  out 
from  Earth  to  rest  at  L2.  Case  4  demonstrated  that  the 
linearized  dynamics  were  valid  for  Case  3.  Therefore  the 
optimal  control  law  successfully  caused  the  spacecraft  to 
be  captured  by  the  moon. 

This  study  only  looked  at  solving  the  two  point 
boundary  value  problem  for  fixed  final  times.  The 
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suggested  next  step  would  be  to  examine  varying  the  final 
time  to  further  optimize  the  trajectory.  This  could  be 
done  by  finding  the  control  Hamiltonian  for  two  different 
final  times  and  then  using  a  secant  method  to  find  the 
final  time  which  renders  the  control  Hamiltonian  equal  to 
zero.  This  might  be  extremely  difficult  since  the  shooting 
method  does  not  converge  for  this  problem  unless  the 
initial  conditions  are  almost  exactly  correct. 
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